%%%%%%%%%%%%%%%%%%%%%%%%
%  analytic solution, also used as boundary function
% 
function [u_1, u_2] = fun_u(xx, yy, E, nu, varargin)
    kappa = 3 - 4*nu;
    G = E/(2 + 2*nu);
    
    r = sqrt(xx.*xx + yy.*yy);  % calculate the radius
    
    % calculate the angle for all points, a little bit tedious!!
    theta = atan(yy./(xx+1e-14));
    idx2 = (xx < 1e-10 & yy > -1e-10);
    theta(idx2) = theta(idx2) + pi;
    
    idx3 = (xx < 1e-10 & yy < 1e-10);
    theta(idx3) = theta(idx3) - pi;
    
    idx = (abs(xx) <= 1e-10);
    theta(idx) = yy(idx)./abs(yy(idx))*pi*0.5;
    
    idx = (r < 1e-10);  % angle for the original (0,0) is
    theta(idx) = 0;     % set to be 0 manually. What a dirty work!
    
    mode = varargin{1};
    switch mode
        case 1   %     % for mode 1:
             lambda = 0.5444837367825;  Q = 0.5430755788367;
             a = kappa - Q*(lambda + 1); 
             b = kappa + Q*(lambda + 1); 
             u_1 = 0.5/G*(r.^lambda).*(a*cos(lambda*theta)  - lambda*cos((lambda - 2)*theta));
             u_2 = 0.5/G*(r.^lambda).*(b*sin(lambda*theta) + lambda*sin((lambda - 2)*theta));
    
    
        case 2    %     % for mode 2:
             lambda = 0.9085291898461;  Q = -0.2189232362488;
             a = kappa - Q*(lambda + 1); 
             b = kappa + Q*(lambda + 1); 
             u_1 =  0.5/G*(r.^lambda).*(a*sin(lambda*theta)  - lambda*sin((lambda - 2)*theta));
             u_2 = -0.5/G*(r.^lambda).*(b*cos(lambda*theta) + lambda*cos((lambda - 2)*theta));

        case 3  %     from paper: Vectorized Matlab codes for lineary Two dimensional 
            alpha = 0.5444837367825;
            omega = 0.75*pi;
            lambda = E*nu/((1+nu)*(1-2*nu)); 
            mu=E/(2*(1+nu));
            
            C1=-cos((alpha+1)*omega)/cos((alpha-1)*omega); 
            C2=2*(lambda + 2*mu)/(lambda + mu); 
            
            ralpha = r.^alpha/(2*mu);
            ur = ralpha.*(-(alpha+1)*cos((alpha+1)*theta) + (C2-alpha-1)*C1*cos((alpha-1)*theta)); 
            ut = ralpha.*( (alpha+1)*sin((alpha+1)*theta) + (C2+alpha-1)*C1*sin((alpha-1)*theta)); 
            
            u_1 = ur.*cos(theta) - ut.*sin(theta);
            u_2 = ur.*sin(theta) + ut.*cos(theta);
       
        otherwise
            error('Do not support Mode %d currently!\n', mode);
    end
    
end